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In this paper, computations of transient, incompressible, turbulent, plane jets using 
the discrete lattice BGK Boltzmann equation are reported. A priori derivation of the 
discrete lattice BGK Boltzmann equation with a spatially and temporally dependent 
relaxation time parameter, which is used to represent the averaged flow field, from its 
corresponding continuous form is given. The averaged behavior of the turbulence field is 
represented by the standard k-e turbulence model and computed using a finite-volume 
scheme on non-uniform grids. Computed results are compared with analytical solutions, 
experimental data and results of other computational methods. Satisfactory agreement 
is shown. 
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1. Introduction 

The computational method based on the lattice Boltzmann equation (LBE) is 
relatively new for fluid dynamics. It is part of the paradigm of simulating complex 
physical phenomena, in particular fluid flows, that are based on the observation 
that the interactions of quasi-particles represented by simple models could give rise 
to very complex emergent phenomena. In 1986, in the seminal works of Frisch, 
Hasslacher and Pomeau^ and Wolfram^, the lattice gas automaton (LGA) was in- 
troduced to simulate the fluid behavior described by the Navier-Stokes equations. 
Their work showed that the key to recovering hydrodynamics from such models is 
that the underlying lattice structure, in which the particles are constrained to move 
and collide while obeying mass and momentum conservation laws, satisfy certain 
symmetry properties. The lattice Boltzmann equation (LBE) was introduced by 
MacNamara and Zanetti'^ to overcome certain drawbacks of the LGA such as the 
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presence of statistical noise and the lack of Galilean invariance. A number of re- 
finements were made, such as the simplification of the description of the collision of 
particle populations^'^ by means of the well-known Bhatnagar-Gross-Krook (BGK) 
model^ which resulted in a considerable simplification of the LBE. In the above 
sense, in which the LBE represents the strcam-and-coUidc picture of the particle 
populations, it may be essentially considered as an extension of the LGA. 

On the other hand, it was believed that the LBE could also be connected to 
the Boltzmann equation, a well known kinetic equation in non-equilibrium statis- 
tical mechanics, that describes the evolution of the particle populations in terms 
of the distribution functions. It was formally shown that the LBE can be derived 
d priori from the Boltzmann equation when its continuous velocity space is con- 
siderably simplified to a certain discrete velocity spacc^^''. With the foundations 
of the kinetic theory, many ideas pertinent to the Boltzmann equation have been 
extended to the LBE or the discrete lattice BGK Boltzmann equation. The rapid 
development and the applications of the lattice Boltzmann method (LBM) that in- 
volves the solution of the LBE has been documented in several review papers^"^^^ 
and monographs^^'^^. In particular, the ability of the method to model physics 
at a smaller scale makes it a potentially promising tool to simulate fiuid fiows in- 
volving interfaces such as multiphase and multicomponcnt flows and other complex 
flows. Algorithmically, the method involves operations that are explicit and local 
and hence naturally amenable for implementation on almost all types of parallel 
computers. Moreover, in the case of incompressible flows, conventional computa- 
tional methods typically require the solution of a Poisson-type equation for the 
computation of pressure field^^'^^ which may be time consuming and only partly 
parallelizable^^; on the other hand, the pressure is always computed locally through 
an equation of state in the LBM. 

Ever since the beginning of the development of the LGA, the precursor to the 
LBE, there were speculations^" about its application to the simulation of turbulence, 
the most dominant form of fiuid flow that occurs in nature as well as engineering 
applications. Because of the inherent limitations of the LGA, some of which were 
discussed earlier, it was noted that its computational requirements would be far 
more demanding than that of conventional macroscopic methods. However, with 
the introduction of the LBE, with its physical as well as computational advantages, 
direct computation of turbulence based on the LBE became feasible. Martinez et 
al.^^ computed decaying turbulence in a shear layer and assessed the accuracy of the 
LBE by comparison of its results with results obtained by employing the spectral 
method. It was shown that the LBE is almost as fast as the spectral method 
on a serial computer and that it may well be more efficient if parallel computing 
strategies are utilized. Some notable work using the LBE aimed at understanding 
the physics of turbulence arc the studies on the enstrophy cascade rangc^^ and 
the energy inverse cascade range^^ in two-dimensional forced turbulence, the study 
of generalized extended self-similarity in three-dimensional, inhomogeneous shear 
turbulence^^ and the work on three-dimensional turbulent channel flow^^. It is 
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interesting to note that there has been a growing interest in the use of kinetic 
theory based approaches involving certain other forms of the Boltzmann equation 
to represent turbulence^^'^^. While such studies, where the LBE or some other forms 
of the Boltzmann equation were used to resolve all the length and time scales, are 
important from a fundamental point of view, for high-Reynolds number flows of 
engineering interest some form of modeling within the LBE framework is desirable. 

Turbulence modeling efforts within the LBE framework are being pursued through 
two approaches: (i) modeling based on the strict kinetic-theoretic formulations; and, 
(ii) modeling based on traditional concepts for which extensive literature already 
exist. The first approach is a more recent one that involves the application of the 
renormalization group (RG) analysis to the simplified form of the Boltzmann equa- 
tion, such as the LBE, to develop a model for turbulcncc^^. It was found that 
this introduces low Knudsen regime closure, a feature that is peculiar to the kinetic 
equation, which the authors believe could potentially offer new physical insights as 
well as alternative mathematical treatments when compared to the Navier-Stokes 
equations. As this approach is still in its early stages, much work remains to be done. 
While the issues related to the first approach continue to be addressed, it is impor- 
tant to develop practical turbulence models for the LBE so that flows of practical 
interest can be computed within the LBE framework. In this respect, traditional, 
statistical averaged procedures to turbulence modeling have been extended to the 
Lgj^28-30 Essentially, in this second approach, the relaxation time parameter that 
appear in the BGK model for the collision term in the LBE is now considered to be 
a spatially and tern,porally dependent variable, instead of a constant. As a result, the 
expression for viscosity that can be obtained from the Chapman-Enskog multiscale 
expansion^^ can be considered to be the sum of the laminar viscosity of the fluid 
and a spatially and temporally dependent eddy viscosity, which can be modeled 
using any statistical averaged approach. An interesting question is whether we can 
derive this discrete lattice BGK Boltzmann equation with spatially and temporally 
dependent relaxation time parameter d priori from its corresponding continuous 
form. In the next section, we address this question. 

It is important to assess the accuracy of the LBE used in conjunction with a 
turbulence model for flows of practical interest. Recently, it has been used to com- 
pute turbulent flow over an airfoil'^^. In this paper, we consider another application, 
namely, the computation of turbulent jets. Fluid jets are encountered in many en- 
gineering applications such as internal combustion engines, gas turbine combustors, 
industrial spray systems and a variety of other situations. Such flows are almost 
exclusively turbulent, with Reynolds number. Re, of the order of 10^ or greater. 
Understanding the transient process of mixing is important in such situations. Di- 
rect simulations of jets would be limited to relatively low Reynolds numbers, of the 
order of few thousands. Hence it becomes imperative to use a model to represent 
turbulence and thereby compute its effect on the resolved flow field. In this work, 
we use the standard k-e turbulence model'^'^ in conjunction with the LBE, such that 
the latter would represent the unsteady, mean flow field behavior, to study transient 
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incompressible plane jets. The time-marching nature of LBM is naturally beneficial 
for this problem that is inherently transient in nature. 

The rest of the paper is organized as follows: In the next section, i.e. Section 2, 
an d priori derivation of the discrete lattice BGK Boltzmann equation, with a 
spatially and temporally dependent relaxation time, from its continuous version is 
provided. The analysis is an extension of that provided in works of He and Luo*^'^. 
In Section 3, the LBE as applied to simulate incompressible flows is discussed. As it 
has been shown that adequate resolution is important to study jets^^, it is necessary 
to employ non-uniform lattice grids within the LBE framework. Hence the LBE 
extension to non-uniform lattice grids is subsequently described in Section 4. This 
is followed by a discussion in Section 5 on the representation of turbulence within 
the LBE framework. The hybrid numerical scheme for the solution of the LBE 
and turbulence model and the computational conditions for the turbulent jet are 
described in Sections 6 and 7 respectively. In Section 8, the computed results are 
compared with analytical solutions, measurements, and published results from other 
computational methods. Finally, the paper closes with summary and conclusions 
in Section 9. 



2. Analysis 

Consider the Boltzmann equation with the BGK form for the collision term, where 
the relaxation time parameter is taken to be a spatially and temporally dependent 
variable 

^+£.Vf = i — (f-r^). (1) 

m ^ ■' x{x,t) ■> ' ^ ' 

Here /(x, is the single-particle distribution function, ^ is the particle veloc- 
ity, A(x, t) is the spatially and temporally dependent relaxation time, and /^^ is the 
local Maxwellian given by 



2RT 



(2) 



where R is the ideal gas constant, D is the dimension of the space, and p, u and T are 
the mean fluid density, velocity and temperature respectively. Limiting ourselves 
to the case of isothermal flows, the fluid density and velocity are obtained as the 
kinetic moments of the distribution function, i.e. 



pu 



j Udi = J U''d^- (4) 



Equation (1) can be formally written in the form of an inhomogeneous ordinary 
differential equation with variable coefflcients 

'^Z I 1 f _ 1 fea 

dt X{^,t)-' A(x,t)-' ' ^' 
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where ^ = ^ + ^ • V is the streaming operator or the time derivative operator 
along the characteristic direction ^. The above equation can be formally integrated 
over a time step of 5t, i.e. 



= p ■^0 



r^* r«t ds 1 



A(x+4s,t+s) . 



+ e 



Jo A(=c+5t',t+t')/(x,^,t), 



(6) 



where s,t' G [0,St]. Assuming that St is small enough and the equilibrium dis- 
tribution fimction, and the relaxation time parameter, A, arc locally smooth 
functions, the following linear approximation by Taylor expansion may be made: 



where 



r« (x + ^t', 4, t + t') = Gi + G2it' + O (S^) , 
A (x + ^t', t + t') = Ai + Aait' + O (S^) , 



(7) 
(8) 



Gi = r (x,^,t),G2i = 



Ai = A (x, t) , A21 



A(x + ^5t,t + 5t)-A(x,i) 



In the above, and in what follows, considering only first order accurate approxima- 
tions in time*"'*, the leading terms of the order of 0{6f) are neglected. Hence, the 
following approximation in the integrand in Eq. (6) may be made: 



1 



A(x + ^i',t + i') + ^ 



Gi + G2it' 1 



A. 



/•*' ds f*' 1 (■,_^2lA^. t' 



Substituting the above in Eq. (6), we get 



it. r^' t!_ 

\l / gAi 

Jo 



Ai Ai 



Gi- 



Ai 



dt' 



+ e-^/(x,^,i). 



(9) 



Now, as 



j e'^dt' = Ai [eT^ ~y'J = ^1 - Ai (e'^ - Ij 
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we can rewrite Eq. (9) as follows: 

Ai f - 1 



1 



= e 



■5* 



Ai 



1 

aT 



A21 
Ai 



Ai { dte'^ - Ai 



+ e-^/(x,^,f). 



Employing Taylor expansion we get e^'^i = 1 ± + O (5^). Substituting this in 
the above equation 



+ 



1 - 
1 - 



A. 

Ai 

A. 

Ai 



-(It 

/(x,4,t), 



1 

aT 



A21 
Ai / V Ai 



and neglecting terms of the order of 0{Sf) 



A 



Ai 



(11) 



(12) 



or, finally we obtain the time-discrete version of the Boltzmann equation 

1 



f{^ + ^6t,i,t + 6t)-f{^,i,t) = - 



(x,t) 



[f{^,^,t)-r{^,u)], (13) 



where T (x, t) = \ (x, t) /St- As shown by He and Luo^'^, the equilibrium distribu- 
tion function may be represented by a truncated small velocity expansion and the 
phase space is discretized such that the numerical quadrature that is used in the 
calculation of the kinetic moments for the conserved variables is exact. Thus, for 
example, in the case of the two-dimensional, nine-velocity (D2Q9) modcl^, shown 
in Fig. 1, the numerical quadrature naturally corresponds to the third-order Gauss- 
Hermite quadrature. Hence, Eqs. (2-4) and (13) may be written in discretized form 
as 



fa {x + eadt,t + 6t) - fa (x, t) = -- 



T(x,t) 

p = ^/„(x,f)=^/^'(x,f), 

a a 

PU = ^fa{^,t)ea = ^f^'' {x,t)ec„ 



[/„(x,f)-/^«(x,t)], (14) 



fS" = WaP 



72 (ea • u) + TTT (ea • u) -—,u 



where 



2c4 

4/9 a = 
«;„=•( 1/9 a= 1,2,3,4 
1/36 a = 5, 6, 7, 8 



2c2 



(15) 
(16) 

(17) 
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and Ba is the discrete set of the velocity space of ^ shown in Fig. 1, with the Cartesian 
component of speed given hy c = 6x/ St where Sx is the lattice spacing and a is the 
velocity direction. Also, fa{x,t) = Wafa{'^,ea,t), /a'^(x,f) = Wafa' {^,ea,t). 
Thus, the discrete lattice BGK Boltzmann equation with spatially and temporally 
dependent relaxation time parameter, similar to the standard discrete lattice BGK 
Boltzmann equation with constant relaxation time parameter, follows d priori from 
its corresponding continuous version and is independent of the LGA. It can be shown 
using the Chapman-Enskog multiscale expansion^^''"'^ that in the long-wavelength 
limit, the mean density and velocity obey the unsteady Reynolds averaged Navier- 
Stokes (RANS) equations with the viscosity related to the lattice parameters, i.e. 

u{x, t) = Uiam + Z^eddy(x, t) = ^C^ ^r(x, t) - ^ St, (18) 

where t'lam is the kinematic viscosity of the fluid and t'eddy(x, t) is the eddy viscosity 
which models the effect of turbulence on the flow field; The calculation of eddy 
viscosity is discussed in Section 5. In addition, it can be shown that the pressure is 
related to density by means of an equation of state 

p(x,i) =c^p(x,t), (19) 

where Cg is the speed of sound, which is equal to c/\/3. 



3. The Incompressible Model 

It may be noted that the LBE as discussed above always simulates the weakly 
compressible RANS equations which is valid for small Mach numbers. Ma. This is 
because of the fact that it models physics locally. Although true incomprcssibility, 
which amounts to infinite speed of sound, cannot be achieved in this model, it 
can be modified such that it minimizes the compressibility eflfects to approximately 
represent incompressible flows. In this work, we employ the "incompressible" LBE 
modcl^''' . According to this model, to approximately represent incompressible flows 
with a constant density po, terms of the order of o(Ma^) in the formulation of the 
LBE are systematically eliminated. This leads to a re-deflnition of the equilibrium 
distribution function, Eq. (17) and a modified expression for the calculation of fiuid 
velocity pou = J2a fa^a- We use the modified equilibrium distribution function as 
derived by He and Luo^^: 



3 9 3 

P + Po <! ^ (ea • u) + ^ (e„ • u)' - 



(20) 



where the coefficient Wa is the same as that used in Eq. (17). In the literature, 
the LBE that results from these modifications is referred to as the "incompressible" 
LBE model. 



4. Non-uniform Lattice Grids 



K.N. Premnath and J. Abraham 



The use of non-uniform lattice grids is desirable in many applications and is impor- 
tant in the simulation of jets where sharp gradients necessitate the use of high reso- 
lution in the near-orifice region^^. The original LBM is restricted to uniform grids, 
in that the minimum streaming distance of the particle populations in one time step 
is exactly equal to the minimum lattice spacing. In other words, the discretization 
of the configuration and particle velocity spaces are coupled. This lockstep advec- 
tion of particle population is a feature inherited from the LGA and is not necessary 
for the; LBE. Since the LBE is actually a simplified form of the Boltzmann equation 
which can be solved without the coupling of the physical and the particle velocity 
spaces, it can be solved on any mesh. Thus, it was proposed by He et al.^^ that the 
collisions still take place on the lattice grids in the manner discussed in the previous 
section; after collisions, the particle populations move according to their velocities 
60,; although the advected distance of the particle populations may not, in general, 
coincide with the mesh spacing, the distribution functions in these locations can 
always be computed using interpolation; after interpolation, the collision and the 
streaming steps are repeated. It has been shown that, if the interpolation method 
is at least of second order, the Navier- Stokes equations can still be recovered from 
the 

In this paper, we employ a second order Lagrange interpolation scheme to im- 
plement this interpolation-supplemented LBE. Figure 2 illustrates the use of this 
method on a stretched non-uniform lattice mesh. For example, for the particle ve- 
locity direction a = 1, at locations which corresponds to the lattice site A, 
{i — 1, j) and {i — 2, j), the collision step is first computed. Then, the particle pop- 
ulations from these locations stream in the positive i direction to a distance |ea|(5t, 
which may not, in general, be equal to the local mesh spacing. After streaming, 
the distribution functions from these three new locations are interpolated to obtain 
the value of the distribution at the location A. For direction a = 3, the lattice 
site locations are taken from the negative i direction as indicated at site B in the 
figure. A similar procedure is adopted for the other particle directions. With this 
arrangement, the incompressible LBE model is used to simulate jets. 

5. Turbulence Modeling 

The standard two-equation k-e turbulence model is employed to represent the effects 
of the length and time scales in the turbulent fiow. It was suggested earlier by Succi 
et al."^^ that these equations might be solved within the LBE framework by creating 
two additional populations, with components in the same directions as the particle 
distribution for each of the two scalar fields. An alternative approach is to solve 
the k-e equations using entirely different computational grids with an appropriate 
numerical scheme'^*'. Here, we consider this latter approach to model the unresolved 
length and time scales in the turbulent, plane jets. The equations for the turbulent 
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kinetic energy, k and the dissipation rate, e are given by 



Po 



-+u.VM = 



d 
dxi 



(IT 



Mlam H 

(Tfc / dXj 



dk 



Tij Sij 



Poe, 



Po 



9a;,- 



, Pt 

Plam H 



de 
dxi 



(21) 



(22) 



The Reynolds stress and strain rate tensors, and Sij respectively, are related by 
the linear constitutive relation through the Bousinessq approximation 



Tij = 2fiTSij - 2/3pokSij, 
and the eddy viscosity is given by 

/ .^ I^T ^ k"^ 

feddy (X, t)=VT = =C^ — 

Po e 



(23) 



(24) 



The values of the model coefficients, Cci,Ce2,o'k,cr£, used in this work are the 
same as that provided in Launder and Spalding^^. The strain rate tensor is directly 
computed from the second kinetic moment of the non-equilibrium part of the dis- 
tribution function, without taking recourse to the finite-differencing of the velocity 
field, using 

^a\^a] {fa fct^^ 



25,r(x,t)4^' 



Po 



(25) 



6. The Hybrid Computational Scheme 

The mean flow field is computed using the LBE, Eq. (14), supplemented by an 
interpolation step, and the turbulence using the k-e model transport equations. 
In this paper, we use the conservative formulation of the k-e equations and solve 
them using a non-uniform finite- volume (FV) based scheme. The computational 
mesh for the k-e equations is the same as that employed for the lattice grids, with 
each finite-volume cell centered at a lattice site, e.g. the location P shown in 
Fig. 3. Following the procedure in Magi^^'^*, the discretized equations are written 
in implicit form with respect to such cells, with convective fiuxes represented by an 
upwind scheme and the diffusive fluxes obtained using central differences across the 
faces, represented by the symbols N , W , S and E in Fig. 3. They are solved using 
the strongly implicit procedure (SIP) due to Stone^^ to obtain rapid convergence of 
the solutions. 



7. Computational Conditions 

Fluid is injected in a plane domain at a constant velocity, JJjnj with ?7inj/c = 0.1, 
where c is the particle speed, through a slot of width d such that d/5x = 4, where 
5x is the minimum lattice spacing, in a relatively large closed chamber. This choice 
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of parameters would keep the flow in the incompressible range so that the applica- 
tion of the incompressible LBE model becomes valid. We consider the density of 
the injected fluid, po, to be equal to 1.0. The ambient fluid, which is also of the 
same density, is assumed to be quiescent initially. The D2Q9 model of Fig. 1 is em- 
ployed. No slip boundary conditions arc imposed at the chamber walls. In this work, 
the boundary conditions are implemented using the non-equilibrium bounce back 
scheme suggested by Zou and He*°. For the k-e equations, the inlet conditions are 
specified, based on assuming equilibrium of turbulence production and dissipation. 
Thus, the inlet turbulent kinetic energy and the dissipation rate are fcin = l.bu'f^ 
and ein = kl^^^/ld respectively. Here, u'^ is the root mean square value of the 
fiuctuating component of velocity and la is the integral length scale. We assume 
u'in = O.Olf/inj and Id = 0.25d. Near the walls, wall functions are used to specify 
the boundary conditions for the turbulent quantities^*. The Reynolds number of 
the jet. Red, is based on the slot width and is defined by Rea = f/injd/fiam- 

Two sets of problems are considered. First, to validate the LBE without a turbu- 
lence model, computations of a laminar plane jet with a Reynolds number of 12 arc 
carried out and compared with similarity solutions. Second, we consider turbulent 
plane jet with Reynolds number of 3 x 10^ and compare the results with measure- 
ments and with those obtained from other computational methods published in the 
literature. We are interested in determining the structure of the jet characterized 
by such quantities as the velocity distribution and the jet half-width, which will be 
defined in the next section. The results will be reported in terms of lattice units, i.e. 
the velocities are scaled by the particle velocity and the distance by the minimum 
lattice spacing. 

8. Results cind Discussion 

When the injected jet travels downstream of the slot, with x as the axial distance 
from the slot and y as the transverse distance from the centerline of the jet, the 

centerlinc axial velocity of the jet progressively decreases as a result of diffusion 
of momentum along the direction transverse to the injected direction of the jet. 
Figure 4 shows the variation of the centerline axial velocity, Uc{x), of the laminar 
jet, normalized by the injection velocity, Umj, as a function of the axial distance 
from the slot, x, normalized by the slot width, d, at different times. It may be 
seen that the velocity distribution is transient in character. Parts of the velocity 
profiles, nevertheless, progressively approach steady state as time progresses. Thus, 
for instance, after 50, 000 time steps, the centerline velocity distribution up to at 
least 60 slot widths downstream of the slot may be considered to be steady. Now, 
the analytical solution for the structure of the jet^^'^^ is applicable only for the 
steady part of the jet and hence they should only be used for comparison with the 
computed results of the steady portion of the jet. The analytical solution based 
on similarity considerations yields the axial velocity, U{x,y), as a function of the 
transverse distance. It is given by the expression U{x, y)/Uc\ = sech rj, where Uci{x) 
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is the centerline velocity and rj is the similarity variable, a dimensionless coordinate 
given by ?7 = 0.5 (Minj/Spo^'iam)^^^ y/a;^/'^. Here, Mi„j is the injected momentum 
flow rate of the jet. The half- width of the jet, 2/1/2, which is defined as the distance 
in the transverse direction from the centerline where U (a;, 2/1/2) = l/2?7ci {x), also 
follows from the similarity solution as yi/2/d = ARe^'^^^ {x/d)^^^ , where A is a 
constant equal to 3.2038. 

Figure 5 shows the variation of the computed half-widths at different times, 
shown in symbols, as a function of the axial distance. Both axes are normalized by 
the slot width. Also plotted on the same figure is the similarity solution. It may be 
seen that when t = 40, 000, at least up to 15 slot widths downstream of the slot the 
computed half-width profile is steady. In this range, the growth of the half-width 
closely agrees with the analytical solution within 4%. 

The computed normalized velocity profile, Uci{x,y)/Uci{x) as a function of 
the similarity velocity variable, 77 at 20 slot widths downstream of the jet when 
t = 60, 000 is compared with the similarity profile in Fig. 6. The computed and an- 
alytical results agree within 3% for r] < 1.5, but the differences increase at distances 
greater than this, due to wall effects. 

Figure 7 shows the lattice/finite- volume mesh employed for the turbulent jet 
computations. A nominal resolution with 400 lattice grids in both the axial and 
transverse directions are employed. The grids are spaced non-uniformly with stretch- 
ing applied in both the directions. The mesh is considered to be uniform in the 
near-field region of the slot, where strong gradients in the flow field is expected. Ef- 
fectively, the mesh accommodates a rectangular domain, whose length in the axial 
distance equals 475d and that in the transverse direction equals 270rf. 

Figure 8 shows the velocity vector field in the domain of the turbulent jet after 
10, 000 time steps. Also shown in the same figure is a plot representing the velocity 
vector field in the near-field of the slot. It may be seen that a pair of vortices, the 
starting vortex or the head vortex pair, are formed from the two shear laj^crs of the 
injected jet and the ambient fluid. As time progresses, it was observed that these 
vortices convect downstream along with the unsteady part of the jet. Also seen 
in the figure is the potential core, in which the jet preserves the injected velocity 
for certain distance downstream from the slot. These basic feature of the jet arc 
consistent with those observed in experiments and in the numerical simulations 
based on the Navier-Stokes equations. 

The computed axial velocity profile, normalized by the centerline velocity, is 
plotted as a function of the transverse distance normalized by the jet half-width at 
20 slot widths downstream of the slot for different times in Fig. 9. It is well known 
that the steady part of the transient jet exhibits asymptotic self-similarity after a 
certain distance downstream of the slot. It may be seen that at the downstream 
location plotted in Fig. 9, the velocity profiles at different times almost overlap with 
one another, implying steady state there. Shown on the same plot in symbols are 
the data from the measurements by Gutmark and Wygnanski^^. The computed 
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results agree within 10% with the measured values. 

As a consequence of self-similarity, the axial variation of the centerline velocity 
and the half-width of the jet obey the following scaling laws^^: Ud (x) ^ x~^l'^ 
and 2/1/2 {x) x. Thus, in the region of self-similarity, (f7inj/?7ci(a;))^ should be a 
linear function of the axial distance, x. Figure 10 shows the computed values of 
(C/inj/t/ci(a;)) , in filled circles, as a function of the normalized axial distance after 
80, 000 time steps. It may be seen that the computed results indeed follow a linear 
variation as expected and agree qualitatively with the measurements. Quantita- 
tively, the data may be fitted to a linear curve that may be represented by 

In this work, this linear curve fit is applied to the computed results for the range of 
distances given by 20 < cc/d < 140 and the constants C\ and Ci thus obtained are 
reported in Table 1, where data from other sources are also compiled for comparison. 
The sources include the three measurements noted above and the direct numerical 
simulation (DNS) data by Stanley et al.'^^ and the large eddy simulation (LES) 
results by Le Ribault et al.'^'^ It may be seen that there is a considerable scatter 
in the values of the constants. In particular, the constant C2 which represents 
the intercept of the linear fit vary significantly, which predominantly reflects the 
quantitative difference between the computed variation and the measurement by 
Gutmark and Wygnanski^^. Physically, this constant is associated with the location 
of the virtual origin of the jet and hence it is expected that different methods could 
result in different values of this constant. 

Figure 11 shows the computed normalized half- width, in filled circles, plotted as 
a function of the normalized axial distance. It may be seen that the computed results 
show a linear trend consistent with the similarity scaling law. For a quantitative 
comparison, the results are fitted according to the following linear curve 

(27) 

The constants Ki and K2 are presented in Table 1. The slope of this curve, i.e. 
dyi/2/dx or Ki is referred to as the spreading rate of the jet and is a constant in the 
similarity region. It is important to reproduce this quantity with sufficient accuracy 
in engineering applications, as it is one of the measures of the rate of mixing of the 
injected and the ambient fluids. In Table 1, along with the computed spreading rate, 
the spreading rate based on the steady Reynolds averaged Navier-Stokes (RANS) 
equations with different two-equation models as reported by Wilcox"^^ and also by 
Magi et aZ.^^ based on the unsteady RANS equations with the standard k-e are 
presented. The measured values for this quantity are in the range 0.100 — 0.110. 
The values based on the direct numerical simulation (DNS) and that based on the 
turbulence models, including the large eddy simulation (LES) turbulence models, 
are in the range 0.090—0.110. The computed value based on the LBE/fc-e FV scheme 



ryi/2 



l/2\ 

d J 
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yields a value for the spreading rate as 0.0965 which is within 5% agreement with 
these other studies. 

9. Summary and Conclusions 

In this paper, it is shown that the LBE with a spatially and temporally depen- 
dent relaxation time paper follows a priori from its continuous version. Compu- 
tations of transient incompressible turbulent plane jets are reported by employing 
this LBE in conjunction with the standard k-e turbulence model equations. The 
turbulence transport equations are solved using a finite-volume (FV) scheme on 
non-uniform grids. The computed structure of the turbulent jet is shown to be con- 
sistent with prior measurements and computed results. In particular, this hybrid 
LBE-FV k-e scheme is found to reproduce the similarity scaling laws for turbu- 
lent jets, i.e., the x~^^'^ decay for the ccntcrlinc jet velocity and linear growth of 
the jet half-width. The computed results agree with the measurements to within 
10%. Laminar jet computations are shown to be in agreement with the analytical 
solution. Although not shown here, we have shown elsewhere^® that this hybrid 
approach is computationally more efficient on parallel computers as compared to 
the conventional schemes^^ when its inherent parallelism is exploited. 
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Table 1. Comparison of the coefficients of the fitted Uneax curve for the normalized inverse square 
axial velocity and the normalized axial variation of the half- width of the turbulent jet. 



Source 


Red 


Ci 


C2 


Ki 


K2 


LBE / Std. k — e model - this work 


30,000 


0.1371 


0.0607 


0.0965 


-0.488 


Unsteady RANS / Std. k - t model** 








0.103 




Steady RANS / Std. k-e modcl^s 








0.108 




Steady RANS / Std. k - w model^^ 








0.101 




LES/dynamic Smagorinsky SGS model*'^ 


3,000 


0.100 


0.89 


0.094 


1.38 


LES/dynamic mixed SGS model^^ 


3,000 


0.220 


0.18 


0.106 


0.40 


DNS^s 


3,000 


0.201 


1.23 


0.092 


2.63 


Measurement - Thomas & Chu (1989)"^ 


8,300 


0.220 


-1.19 


0.110 


0.14 


Measurement - Browne at o/.(1983)** 


7,620 


0.143 


-9.00 


0.104 


-5.00 


Measurement - Gutmark & Wygnanski(1976)*^ 


30,000 


0.123 


4.47 


0.100 


-2.00 



K.N. Premnath and J. Abraham 




Fig. 1. D2Q9 Lattice. 



Turbulence modeling of jets in discrete lattice BGK Boltzmann method 



1 1 1 1 i 1 1 1 1 i- 



->- 



B 



j + 2 



j+1 



j-2 



1 1 t t t t t t — t — i 



i -2 i -1 



i i+1 



i +2 



> I 



Fig. 2. Non-uniform lattice mesh. 
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Fig. 3. Finite control volumes embedded on lattice mesh. 
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Fig. 4. Transient axial velocity profiles as a function of axial distance of the laminar jet; Red = 12. 
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Fig. 5. Comparison of computed results and similarity solution for half-width of the laminar jet; 
Rea = 12. 
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Fig. 7. Computational lattice/finite volume mesh employed. 
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Fig. 8. (a) Velocity vector field in the domain of the turbulent jet after 10, 000 time steps; (b) 
Velocity vector field in the near-field of the slot showing the potential core and the starting vortex 
pair; Re^ = 30 X 10"*. 
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Fig. 9. Comparison of computed results and measurements for the normalized axial velocity profile 
as a function of the similarity variable of the turbulent jet; iie^ = 30 x 10**, x/d = 20. 
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Fig. 10. Normalized inverse square of the centerline velocity as a function of the normalized axial 
distance of the turbulent jet; Re^ = 30 x 10^. 
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Fig. 11. Normalized jet half-width as a function of the normalized axial distance of the turbulent 
jet; Red = 30 x 10*. 



